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ABSTRACT 

We present a set of new numerical methods that are relevant to calculating radiation 
pressure terms in hydrodynamics calculations, with a particular focus on massive star 
formation. The radiation force is determined from a Monte Carlo estimator and enables 
a complete treatment of the detailed microphysics, including polychromatic radiation 
and anisotropic scattering, in both the free-streaming and optically-thick limits. Since 
the new method is computationally demanding we have developed two new methods 
that speed up the algorithm. The first is a photon packet splitting algorithm that 
enables efficient treatment of the Monte Carlo process in very optically thick regions. 
The second is a parallelisation method that distributes the Monte Carlo workload 
over many instances of the hydrodynamic domain, resulting in excellent scaling of the 
radiation step. We also describe the implementation of a sink particle method that 
enables us to follow the accretion onto, and the growth of, the protostars. We detail 
the results of extensive testing and benchmarking of the new algorithms. 
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1 INTRODUCTION 

Massive stars are hugely important in galactic ecology, en¬ 
riching them chemically and providing strong feedback ef¬ 
fects via radiation, stellar winds and supernovae. They also 
play a pivotal role in measuring star formation rates in dis¬ 
tant galaxies, and therefore determining the star formation 
history of the Universe ( Kennicutt [1998 1. However, the pro¬ 
cess by which massive stars form is less well understood than 
that of lower-mass stars, both observationally, because mas¬ 
sive protostars are rare and difficult to observe, and theoret¬ 
ically because radiation feedback has a much stronger influ¬ 
ence on the gas dynamics than it does for solar-mass objects. 
Broadly it appears that observational and theoretical per¬ 
spectives on massive star formation are converging towards 
a scenario in which a massive young stellar object contains a 
central (> 10 M©) protostar surrounded by a stochastically 
accreting Keplerian disc and a stellar- or disc-wind outflow, 
although substantial conflicts between the predictions of nu¬ 
merical models remain | Tan et al.|2014 |. 

It has been known for some time that massive proto¬ 
stars undergoing spherical accretion will eventually produce 
sufficient radiation pressure to halt the inflow of material 
(e.g. Wolfire fc Cassinelli||1987 1. Extending models to mul¬ 


tiple dimensions, and the introduction of angular momen¬ 
tum, leads to a scenario in which the protostar can accrete 
via a disc, which intercepts a smaller fraction of the pro- 
tostellar radiation field, and the radiation escapes preferen¬ 
tially via low-density cavities along the rotation axis (e.g. 
Yorke fc Sonnhalter|2002 1. Recent work by |Krumholz et ^ 


120091 used 3-D adaptive mesh rehnement (AMR) simula¬ 
tions coupled with a flux-limited diffusion (FLD) approxima¬ 
tion for the radiation transport (RT) to simulate the forma¬ 
tion of a massive binary. They found accretion occurred via 
a radiatively-driven Rayleigh-Taylor instability that allowed 
material to penetrate the radiation-driven cavities above and 
below the protostars. Their simulation resulted in the for¬ 
mation of a massive (41M0and 29 M©) binary pair. 


Kuiper et al. (20101 implemented a hybrid RT scheme 
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in which ray tracing is used to calculate the radiation field 
from the protostellar photosphere out to where the material 
becomes optically thick, at which point the ray-traced field 
becomes a source term in the grey diffusion approximation 
used for the rest of the domain. The method shows good 
agreement with more detailed treatments of the RT but has 
vastly reduced computational cost ( [Kuiper k, Klessen[2013P . 
Two-dimensional calculations adopting this method resulted 
in the formation of massive stars via stochastic disc accretion 
accompanied by strong radiatively driven outflows ( [Kuiperj 
et al.|2dl0 1. 
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2 Tim J. Harries 


Radiatively-driven Rayleigh-Taylor instabilities have 


not been identified in the hybrid RT models (Kuiper 


et al. 20121, suggesting a fundamental difference with the 


Krumholz et al. computations. It appears that the grey FLD 
method considerably underestimates the driving force of the 
radiation, since the opacity used for the radiation pressure 
calculation is the Rosseland opacity at the dust tempera¬ 
ture, whereas in the hybrid case the primary momentum 
deposition occurs where the radiation from the protostar is 
incident on the dust (near the dust sublimation radius). At 
this point the radiation temperature is much higher, and the 
dust is much more opaque. It is possible that the grey FLD 
method underestimates the the driving force by up to two 
orders of magnitude ( Kuiper et al.|[2012 K 

The fundamental point here is that the dynamics of 
the simulations can be critically affected by the level of mi¬ 
crophysical detail employed. The lessons of recent numeri¬ 
cal simulations are that disc accretion may occur, and fast, 
wide-angle outflows may be driven. It seems that radiation 
pressure does not present a fundamental physical barrier to 
massive star formation via accretion, although the strength, 
geometry, and longevity of the outflows is still contentious. 
Furthermore, neither the Krumholz et al. (20091 models 
nor those of the Kuiper et al. (20101 include the effects of 
photoionization, which will come into play as the protostar 
shrinks towards the main sequence and its radiation field 
hardens. All these factors point towards the necessity for 
a more thorough treatment of the microphysical processes 
that underpin that radiation feedback, such as a polychro¬ 
matic prescription for the radiation field, consideration of 
both dust and gas opacity (absorption and scattering), and 
the inclusion of photoionization. 

Recently significant progress has been made in 
using MC transport in RHD codes. For example 
|Nayakshin et al.| ( |2009| ) used MC photon packets 
to treat radiation pressure in smoothed-particle hy¬ 
drodynamics. Acreman et al. (2010) presented a 
proof-of-concept calculation involving combining 3D 
smooth-particle hydrodynamics with MC radiative 
equilibrium. In this case the SPH particle distribu¬ 
tion was used to construct an AMR grid on which 
the MC radiative-equilibrium calculation was con¬ 
ducted. We then developed an Eulerian hydrody¬ 
namics module to perform 3D hydrodynamics on 
the native AMR grid used for the radiation transfer 
(Haworth &: Harries 2012 |. A similar approach was 
adopted by Noebauer et al. (2012), who incorpo¬ 
rated a radiation-pressure force estimator as well as 
a radiative-equilibrium calculation. [Roth &: Kasen] 
(2014) also coupled MC transport and ID hydro¬ 
dynamics, incorporating special relativity and reso¬ 
nance line scattering. 

We have now developed an AMR RHD code (torus) 
that employs a highly-parallelised, time-dependent, Monte- 
Carlo (MC) method ( |Harries||2011[ ) to follow the RT at a 
level of microphysical detail comparable to that of dedicated 
RT codes such as Cloudy ( Ferland et al.|[2013 l. The conse¬ 
quences of this breakthrough are two-fold: firstly, the in¬ 
creased sophistication of the treatment of the radiation field 


can have a significant impact on the hydrodynamics (Ha- 
|worth fc Harries||2012[ ), and secondly it becomes possible 
to make a much more direct comparison between the mod¬ 


els and observations across a wide range of continuum (e.g. 
near-IR dust, radio free-free) and line (e.g. forbidden, molec¬ 
ular or recombination) diagnostics, ( Rundle et al.|20T0 Ha- 
worth et al.|201^ [2013 (. 


In this paper we describe the implementation of (i) a 
new, highly parallelised method for calculating the radiation 
pressure using Monte Carlo estimators, (ii) a packet splitting 
method that improves the efficiency of the method at high 
optical depths, and (iii) a Lagrangian sink-particle algorithm 
that enables us to follow the growth of the protostar. We 
provide the results of extensive benchmarking of the new 
methods, with a view to conducting simulations of massive 
star formation that incorporate both radiation pressure and 
ionisation feedback. 


2 RADIATION PRESSURE 


Our Monte Carlo method is closely based on that of |Lucy| 
(1999), in which the luminosity L of the illuminating ob¬ 
ject is divided up into N discrete photon packets that have 
a constant energy as they propagate through the adaptive 
mesh. The energy of each packet is 


= WiLAt 


( 1 ) 


where At is the duration of the Monte Carlo experiment and 
Wi is a weighting factor, normalized so that 

N 

= ( 2 ) 

i = l 


Note that the original Lucy method has Wi = 1/N, so that 
each photon packet had the same energy. 

As the packets propagate through the grid they may un¬ 
dergo absorption or scattering events. During an absorption 
event the photon packet is immediately re-emitted, with a 
new frequency selected at random from a probability density 
function created from the appropriate emissivity spectrum 
at that point in the grid. The radiation field therefore re¬ 
mains divergence free during the calculation. After the N 
packets have been propagated, the energy density for each 
cell can be straightforwardly computed and hence the ab¬ 
sorption rate can be determined. 

The calculation of the radiation pressure may be conve¬ 
niently conducted in parallel with the radiative equilibrium 
calculation. The simplest method is to assign a momentum 
change to each cell (Apj), which is zeroed at the start of 
each MC loop. The each time a photon packet enters and 
then leaves a cell the difference in momentum between the 
incoming and outgoing momenta is added to the cell’s mo¬ 
mentum change 


Apj 


Apj -f 


e* 

c 


(Uin Uout) 


( 3 ) 


where Uin is the direction vector of the photon packet when it 
enters the cell and Uout is the direction vector of the photon 
packet when it exits the cell. At the end of the MC loop we 
can then calculate the radiation-force per unit volume for 
each cell via 


frad.j 


Apj 

AtV, 


( 4 ) 


where V) is the cell volume. We will call this the momentum 
tracking algorithm. 
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Massive star formation simulation algorithms 3 


This method has the advantage that it is straightfor¬ 
ward to implement, is transparently related to the underly¬ 
ing physics, is valid for arbitrary scattering phase matrices, 
and is very fast to calculate. It does however have the disad¬ 
vantage that the number of scatterings/absorptions per cell 
tends to zero in the optically thin limit. This is an analo¬ 
gous problem to that associated with the |Bjorkman fc Wood| 
(20011 radiative-equilibrium method. 

An alternative method is to use a MC estimator of the 
radiation vector flux to calculate the radiation force. We 
define the energy density Ui, as the energy density per unit 
frequency of a beam of radiation. The energy carried in a 
beam of volume dV is 


dE = UudV dv (5) 

and since dV — cdAdt we get 

dE = Ui^cdAdtdu = EdtdudAdQ, (6) 


where E is the specific intensity. Hence we have 

w, = —dn. (7) 

c 

If a photon packet traverses length £ between successive 
events (events being absorptions/scatterings or crossing cell 
boundaries) then the energy density due to this part of the 
path of packet Si is {ei/Vj)St/At where St = Ijc. We can 
therefore obtain an MC estimator of the energy density 


Uvdu 


ed 

VjcAt 


( 8 ) 


This may be related to the specific intensity via equation 
and we get 


IudQ.dv = 


ed 


( 9 ) 


The radiation flux vector of cell j may then be written 


= ( 10 ) 

du 


where u is the unit vector in the direction of the photon 
packet as it traverses £ and only the paths for photon packets 
with frequencies in the range {v, u + dv) are considered. The 
force per unit volume for cell j is then 


fradjii 


K^pFiydv = 


1 1 
c At 




( 11 ) 


where is the opacity at the frequency v of packet Ci as it 
traverses length £, p is the mass density, and the summation 
is now over all packets. Note that the opacity must consist 
of both the absorption and scattering opacities 


Ki = {1 — a^)Ki, -I- (1 — (12) 

where is the albedo, and = {cos 9) where 6 is the 
scattering angle (i.e. gi, = 0 for isotropic scattering, positive 
for preferentially forward scattering and negative for back 
scattering). Although marginally more expensive to calcu¬ 
late, this flux estimator method should be superior in the 
optically-thin limit, since even photon packets that do not 
interact in a cell contribute to the estimator. 


2.1 Radiation pressure tests 

We constructed a suite of tests of our radiation pressure 
algorithms based on a uniform-density sphere of radius 
Rs = 0.1 pc and mass 100 M© illuminated by a solar- 
type star. The grey opacity k was selected in order that 
r = KpRs = 0.1, 1,10 or 100. The models were run with a 
1-D spherical geometry on a uniform radial mesh comprising 
1024 cells. We ran models for pure absorption (a = 0) and 
pure scattering (a = 1) cases, and for the scattering models 
we further subdivided the models into isotropic scattering 
and forward scattering (using a Henyey-Greenstein phase 
function with g — 0.9). We used 10® photon packets for 
all the models. For each model we calculated the radiation 
force per unit volume (Fnp/c) for both the Monte Carlo flux 
estimator and the momentum tracking algorithm. 


2.1.1 Pure absorption case 

In this case the radiation force per unit volume (/rad) as a 
function of radius r in the sphere is expected to be 

LKpexp(-T(r)) 

where r(r) is the radial optical depth from the centre of 
the grid to r. Since the photon packets stop propagating 
when they are absorbed the number of packets passing each 
cell will be monotonically declining as a function of radius, 
and we therefore expect the noise on the estimates of the 
radiation force to increase radially. Furthermore we expect 
that the noise in the momentum tracking algorithm to be 
larger than that of the flux estimator method, particularly 
for the optically thin cases (r = 0.1,1). 

The results of this benchmark are displayed in Figure[^ 
and the Monte Carlo estimators show good agreement with 
the analytical result, with the expected noise dependencies. 
These models ran rapidly, and since each photon packet only 
undergoes one absorption event the higher optical depth 
cases ran the most quickly. 


2.1.2 Pure scattering cases 

For isotropic scattering the number of photon interactions 
before escape from an optically-thick, uniform spherical of 
radial optical depth Tmax will be 

A'acat « Tmax/(T^) = T-iax/2. (14) 


Thus for the r = 100 model each photon packet will undergo 
~ 5000 scattering events before escaping. Since the radiation 
field is divergence free, the outward flux at every radius is 
a constant F — LKinr^) and the radiation force per unit 
volume is then 


_ Lkp 
47rr2c' 


(15) 


Once again we expect the flux estimator method to have 
less noise than the momentum tracking method, and that 
the latter algorithm will be rather poorer in the low optical 
depth models. 

The results of the isotropic, pure scattering benchmark 
are displayed in Figure Excellent agreement is seen be¬ 
tween the Monte Carlo estimators and the analytical result. 
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4 Tim J. Harries 




Figure 1. Benchmark for absorption {a = 0). The solid lines 
indicate the analytical results (FKp/c), while values calculated 
from the Monte Carlo flux estimator are shown as plusses (+) 
and those found from following packet momenta are shown as 
crosses (x). Benchmarks are plotted for spheres of four optical 
depths (r = 0.1,1,10,100). 



R/Rs 


Figure 2. Benchmarks for an isotropic scattering medium (a = 
1). Symbols are the same as those for Figure^ 


and as expected the noisiest estimator is the momentum 
tracking algorithm for the r = 0.1 model. 

Finally we ran the pure scattering model for a forward¬ 
scattering Henyey-Greenstein phase function (g = 0.9). In 
this model the nett radiation force should be reduced by a 
factor of (1 — 5 ') over the corresponding isotropic scattering 
case, i.e. 


Figure 3. Benchmarks for a scattering medium (a = 1) with a 
forward-scattering Henyey-Greenstein phase function {g = 0.9). 
Symbols are the same as those for Figure^ 


Lnpjl - g) 


(16) 


The results of the forward-scattering model are given in Fig¬ 
ure and excellent agreement with the expected analytical 
radiation force is found. 


2.2 Radiation pressure hydrodynamic test 


Satisfied that the radiation pressure forces were being prop¬ 
erly captured by the code, we progressed to testing the treat¬ 
ment in a dynamic model using a similar benchmark to that 
presented by Nayakshin et al. (20091. A uniform density 
sphere (po = 1-67 x 10“ gcc“ and Rs = 10^^cm) was 
modelled using 1-d grid comprising 1024 uniformly spaced 
radial cells. The sphere was illuminated by a 1 Lq star placed 
at the centre, and we used 10® photon packets per radiation 
step. 


For the purposes of testing we assumed isothermal gas 
at a temperature of 10 K and considered cells above a critical 
density threshold (1.67 x 10“^® gcc“^) to be completely op¬ 
tically thick at all wavelengths with an albedo of zero, while 
cells with densities below the threshold were considered to 
be transparent to all radiation. Photon packets entering op¬ 
tically thick cells were immediately absorbed and were not 
remitted, allowing straightforward comparison with anal 3 rt- 
ical results. 

In this case the radiation pressure sweeps up a thin shell 
whose equation of motion (assuming a negligible contribu¬ 
tion from thermal gas pressure) is given by 



(17) 


Integrating the above expression we get 


R = 


3L 

2ncpo 


1/4 


(18) 
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Radius (pc) 

Figure 4. The development of the radiation-pressure driven 
spherical shell as a function of time. The mass density as a func¬ 
tion of radius is plotted at 1000 kyr intervals from a simulation 
time of 1000 kyr to 8000 kyr. Some broadening of the shell is evi¬ 
dent at later times. 



Figure 5. The evolution of the shell radius with time for the 
radiation pressure hydrodynamic test (solid line). The theoreti¬ 
cal shell expansion based on the thin-shell approximation is also 
shown (dashed line). 


We followed the development of the shell for 8 kyr, writ¬ 
ing out the radial density profile at 10^ s intervals (see Fig¬ 
ure 0, and we then determined the position of the shell 
for each radial profile. This was done by making a parabolic 
ht to the peak density in the grid and the density at the 
adjacent radial grid points. We show the results of the dy¬ 
namical calculation in Figure There is good agreement 
between the model calculation and the theoretical growth 
of the shell predicted by equation Deviations from the¬ 
ory at late times may be attributed to departures from the 
thin shell approximation. 


3 PACKET SPLITTING 


The splitting of particles in the MC method has a 
long history, dating back to the first neutron trans¬ 
port codes ( Cashwell fc Everett||l959 )■ Splitting, and 
its inverse, the so-called Russian Roulette method, 
are variance reduction techniques designed to im¬ 
prove the efficiency with which the MC sampling 


operates. Although the original Lucy ( 19991 algorithm had 


equal energy photon packets, this is not a fundamental re¬ 
striction. In the case of star formation calculations the pho¬ 
ton packets may be produced within gas with very high 
optical depth, and the packets may undergo thousands of 
scattering events prior to emerging from the computational 
domain. The MC estimators of the energy density and ra¬ 
diation pressure will be good quality for the optically thick 
cells where packets spend most of their time, and thus it 
is inefficient to use equal energy packets. Instead we em¬ 
ploy a packet splitting algorithm, in which a lower number 
of higher energy packets (Ahigh) are released from the pro¬ 
tostar. These propagate through the optically thick region 
and then are each split into a number of lower energy packets 
(Alow) in the optically thin region (see|^. The total number 


of packets that emerge from the grid is then N = AhighAiow 
The key is to correctly identify the point at which the packet 
splitting takes place, in order that (i) high energy packets 
do not propagate into optically thin regions, since they will 
increase the noise in the MC estimators, (ii) low energy pack¬ 
ets undergo a small (but non-zero) number of interactions 
before leaving the computational domain. 

We constructed a one-dimensional test model for the 
packet splitting algorithm, comprising a 0.1 pc radius sphere 
containing 100 Mq of gas and an r~^ density profile. The 
dust opacity was from Draine & Lee ( 19841 silicates with a 
dust-to-gas density ratio of 1 per cent. The dust is heated 
by a central luminous protostar of Teff = 4000 K and radius 
R — 150 Rq. We first calculated the radiative equilibrium 
without packet splitting using N = 10® photon packets, and 
used the final iteration of the radiative equilibrium calcula¬ 
tion as our benchmark. We subsequently ran the same model 
with packet splitting and Ahigh = 1000 and Aiow = 100, 
defining the high optical depth region as that for which 
the Planck-mean optical depth to the sphere radius was 
20. Note that in two or three-dimensional problems 
the integral of the Planck-mean optical depth is cal¬ 
culated in the positive and negative directions of 
each coordinate axis, defining a typically ellipsoidal 
boundary for packet splitting. 

We plot the temperature profile of the sphere in Fig- 
ure[^for both cases, and find good agreement. In particular 
there is a smooth change in the temperature through the op¬ 
tical depth boundary (the region where the packet splitting 
occurs). 

In the calculation without splitting each photon packet 
undergoes ~ 20 000 absorption or scattering events prior to 
its emergence from the computational domain. When packet 
splitting is enabled each high energy photon packet still un¬ 
dergoes ~ 20 000 absorptions or scatterings, but the low 


© 2015 RAS, MNRAS 000,[iHTT] 




























6 Tim J. Harries 



Figure 6. A cartoon illustrating the packet splitting method. A 
high energy packet is released from the protostar and undergoes 
many absorption and scattering events (black arrows). Eventually 
the packet passes beyond the region identified as being optically 
thick (dashed line) and at this point the packet is split into A^iow 
low-energy packets (red arrows) which eventually emerge from 
the computational domain (here Ni^-^ = 5). 



R/Rs 


Figure 7. The results of the packet splitting test. The temper¬ 
ature of the sphere is plotted as a function of radial distance 
(in units of the sphere radius Rs = 0.1 pc). The model without 
packet splitting (-1- symbol) shows excellent agreement with the 
temperature found using packet splitting (x symbol). The verti¬ 
cal line indicates the radius beyond which the high energy packets 
undergo splitting. 


energy photon packets only undergo ~ 200. This results in 
a speed-up of the radiation step of a factor of ~ 50. 


4 PARALLELISATION 


The primary drawback in such a detailed treatment of the 
radiation field is the computational effort involved. Fortu¬ 
nately the MC method, in which the photon packets are es¬ 
sentially independent events, may be straightforwardly and 
effectively parallelised. The top level of parallelisation in¬ 
volves domain decomposing the octree that stores the AMR 
mesh. Each sub-domain belongs to a separate MPI thread, 
enabling the code to run on distributed memory machines 
and enabling the use of domains with a larger memory foot¬ 
print than that available on a single node. We choose a sim¬ 
ple domain decomposition, which although not necessarily 
load balancing, does allow a straightforward implementation 
in the code. For eight-way decomposition each branch of the 
octree is stored on an individual MPI process (with an ad¬ 
ditional thread that performs tasks such as passing photon 
packets to threads). Similarly 64-way decomposition may be 
achieved by domain decomposing further down the branches 
of the octree, and this is the decomposition that is employed 
for majority of our runs (although 512-way decomposition 
is implemented we do not have access to the resources nec¬ 
essary to regularly run the code in this mode). 

The main bottleneck is the communication overhead 
when passing photon packets between threads, and we opti¬ 
mize this by passing stacks of photon packet data between 
the MPI threads rather than individual packets in order to 
reduce latency. Thus when a photon packet reaches 
a boundary between domains (naturally this always 
corresponds to cell bonndary) then the packet posi¬ 
tion, direction, frequency, and energy of the photon 
packet is stored on a stack. Once the stack reaches a 
set number of packets (typically 200), then the stack 
is passed to the appropriate MPI thread. We note 
that this algorithm closely resembles the MILAGRO 
algorithm described by Brunner et al. (2006|. Fur¬ 
ther optimisation of the stack size is possible, includ¬ 
ing varying the stack size across the domain bound¬ 
aries and dynamically altering the stack size as the 
calcnlation progresses, but we have yet to implement 
this. 


A further level of parallelisation is achieved by having 
many versions of identical computational domains (which 
we refer to as hydrodynamic sets), over which the photon 
packet loop is split, with the results derived from the radia¬ 
tion calculation (radiation momenta, cell integrals etc) col¬ 
lated and returned to all the sets at the end of each iteration. 
The thermal balance and ionization equilibrium calculations 
(which can be time consuming) are then parallelised over the 
sets, with each thread corresponding to a particular domain 
performing the equilibrium over an appropriate fraction of 
its cells before the results are collated and distributed (see 

We performed a scaling test to demonstrate the effi¬ 
cacy of our scheme. We ran the model detailed in section 
in 2D. The benchmark used one hydrodynamics set with 
the quadtree domain-decomposed over 4 threads, giving a 
total of 5 MPI threads including the control thread. We 
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Figure 8. A schematic showing the parallelisation of the radia¬ 
tion loop in TORUS. Here the photon packets are split across two 
sets of the AMR mesh (green squares) and each mesh is divided 
into four domains (blue squares) , with inter-domain communica¬ 
tion of photon packets (indicated by the black arrows). At the end 
of the iteration the Monte-Carlo estimators are combined for each 
domain in order to compute revised temperatures and ionization 
fractions of the gas. 


Number of MPI threads 


5 10 20 50 



Number of sets 


Figure 9. Results of the scaling test of the parallelisation. The 
same 2-D model was run with 1, 3, 6 and 12 hydrodynamics sets, 
corresponding to 5, 15, 30 and 60 MPI threads. The wall time 
for one radiation timestep is shown (normalized by the 5-thread 
calculation). The solid line indicates perfect scaling. 


recorded the wall time for one radiative transfer step, and 
subsequently ran this same model for 3, 6, and 12 hydro¬ 
dynamics sets (corresponding to 15, 30 and 60 MPI threads 
respectively). We found excellent scaling (see Figure § up 
to 60 threads, with the slight deviation from embarrassing 
scalability due to the all-to-all MPI communication neces¬ 
sary when the results of the MC estimators are gathered 
(this is approximately 10 per cent of the wall time for the 
60 thread run) 


5 SINK PARTICLE IMPLEMENTATION 


We implemented a Lagrangian sink particle scheme within 
TORUS in order to follow the formation of stars in our hy¬ 
drodynamics code. Originally adopted for smoothed-particle 
hydrodynamics solvers ( Bate et al.||1995 ), sink particles al¬ 
low one to overcome the restrictively small Courant time 
associated with high densities as the collapse proceeds, by 
removing gas from the computation and adding it to point¬ 
like particles that interact with the gas via gravity (and 
possibly radiation), but not via thermal gas pressure. 

The fundamental principles for creating a reasonable 
sink particle implementation are (i) creating sink particles 
only as a ‘last resort’, (ii) accreting gas in a manner that has 
minimal impact on the dynamics of the gas immediately 
around the sink particle, (iii) properly accounting for the 
gravitational forces between Lagrangian sink particles and 
the gas on the grid. 

The use of sink particles in grid-based hydrodynamics 
codes has been investigated by [Krumholz et ah] ( |2004[ ) and 
Federrath et al. (20101. Broadly speaking the two implemen¬ 


tations described deal with items (i) and (iii) in the same 
manner, but differ in the way that accretion is treated. The 
Federrath et al. method defines an accretion radius (as a 
small multiple of the smallest cell size in their AMR mesh) 
and simply tests that gas in cells within the accretion radius 
is bound to the sink particle, and then accretes the gas mass 
(and its associated linear and angular momentum) above a 
threshold density onto the sink particle. The Krumholz et al. 
method uses a dynamically varying accretion radius, which 
ranges from 4-cells to the Bondi-Hoyle radius. The method 
then ascribes an accretion rate onto the sink based on Bondi- 
Hoyle-Lyttleton accretion and removes the appropriate mass 
from each cell within the accretion radius based on a weight¬ 
ing function that falls rapidly with distance from the sink 
(cells outside the accretion radius have a weighting of zero). 
The advantage of this method is that an appropriate accre¬ 
tion rate is maintained in situations where the Bondi-Hoyle 
radius is unresolved, whereas the Federrath method breaks 
down in this regime. However the Krumholz method relies 
on a statistical smoothing of the flow in order to determine 
the appropriate far-held density and velocity to use when 
calculating the Bondi-Hoyle accretion rate. In contrast the 
Federrath method relies primarily on the accretion flow be¬ 
ing supersonic outside the accretion radius, and therefore 
the algorithm by which material is removed from the grid 
has no impact on the upstream dynamics. 

In our implementation we use a method akin to Feder- 
rath’s, since we are always in the regime where the Bondi- 
Hoyle radius is well resolved. Below we describe the details 
of our sink particle implementation, and several accretion 
and dynamics tests of the method. 


© 2015 RAS, MNRAS 000,[iHTT] 











































































































8 Tim J. Harries 


5.1 Sink particle creation 


We adopt the same criteria for sink particle creation as|Fed-| 


errath et al. (20101. For the cell under consideration we de¬ 


fine a control volume that contains all cells within a prede¬ 
fined radius race- Before a sink particle is created a number 
of checks on the hydrodynamical state of the gas in the con¬ 
trol volume must be passed, briefly: 


• The central cell of the control volume must have the 
highest level of AMR refinement. 

• The density of the central cell must exceed a predefined 
threshold density pthresh, thus ensuring that V • (pv) < 0 for 
that cell. 

• Flows in cells along the principle axes must be directed 
towards the central cell. 

• The gravitation potential of the central cell must be the 
minimum of all the cells in the control volume. 

• The control volume must be Jeans unstable i.e. 

IFlgravl > 

• The gas must be in a bound state i.e. iSgrav + J?th + 
IFkin < 0, where JFkin is the kinetic energy of the gas in the 
control volume where the speeds are measured relative to 
the velocity of the centre of mass of the control volume. 

• The control volume must not overlap with the accretion 
radius of any pre-existing sink particles. 


If these tests are passed then a sink particle is created at 
the centre of the control volume, and accretes gas according 
to the method detailed below. 


5.2 Accretion onto sink particle 

Each sink particle has an associated accretion radius race 
which is defined as 2.5 times the size of the smallest cell 
in the AMR mesh. The accretion radius defines a spherical 
volume, and after each hydrodynamics step each the cells 
in the accretion volume are checked against a predefined 
density threshold (pthresh). Cells with a density above pthresh 
undergo a further check to see if the gas in them is bound 
to the sink particle, that is 

Egrav J- i?kin + Eth < 0. (19) 

Note the kinetic energy of the gas is measured relative to 
the velocity of the sink particle. If the gas is both bound 
and above the density threshold then the mass of the sink 
particle is increased by (p — pthresh)F where p and V are 
the density and volume of the cell. The linear and angular 
momentum of the accreted mass is added to the sink particle 
and subtracted from the cell. 


5.3 Sink particle motion 

The gravitational influence of the gas on a given sink particle 
is found by summing up the gravitational forces from all the 
cells in the AMR grid. For all cells excepting that containing 
the sink particle we use 

Fj = ^ GMjpiVis{rij,h)i (20) 

i 


and the sink particle respectively, h is the gravitational soft¬ 
ening length, and sifij, h) is a cubic spine softening function 
(see equation Al of Price fc Monaghan||2007 (. 

For the cell containing the sink particle we instead per¬ 
form a sub-grid calculation by splitting the mass of the cell 
into 8 ® subcells and sum the force over the subcells in an 
analogous manner equation | 20 [ 

The sink-sink interaction is computed using a sum of the 
gravitational forces over all sinks (this is computationally 
tractable since the number of sink particles is small). The 
equation of motion of the sink particles is integrated over a 


timestep by using the Bulirsch-Stoer method (Press et al. 
T9^, once again using the cubic spline softening of|Price| 


& Monaghan (20071. 


5.4 Gas motion 

The self-gravity of the gas is solved using a multigrid solu¬ 
tion to Poisson’s equation, and with the gradient of the po¬ 
tential appearing as a source term in hydrodynamics equa¬ 
tions. 


6 SINK PARTICLE DYNAMICS TESTS 


Following Rubber et al. (20111 we adopt the three-body 
Pythagorean test problem of Burrau| (19131 to benchmark 
the sink-sink gravitational interactions and test the integra¬ 
tor. This problem has masses (taking G = 1) of mi = 3, 
m 2 = 4 and m 3 = 5 starting at rest at Cartesian coordi¬ 
nates (1,3), (-2,-1) and (1,-1) respectively. The subse¬ 
quent motion of the masses involves a number of close inter¬ 
actions between sinks, and was first numerically integrated 


by Szebehely & Peters (19671. The results of the test prob¬ 
lem are given in Figure m and show excellent agreement 
with Figures 2, 3, and 4 of Szebehely & Peters (19671. The 


total energy of the system over duration of the computation 
is conserved to better than 1 part in 10 ®. 

The implementation of the gas-on-sink gravitational 
force calculation was tested by using a power-law (p oc 
p{r)~^) density sphere of 100 Mq and radius 0.1 pc, on a 
3D AMR mesh with minimum cell depth 6 (equivalent to a 
fixed-grid resolution of 64 x 64 x 64) and maximum cell depth 
10 (1024 X 1024 X 1024). The computational domain was di¬ 
vided over 64 MPI threads. Three test sink particles of neg¬ 
ligible mass were placed in the grid at radii of 2.5 x 10^^ cm, 
1 X 10^^ cm, and 5 x 10^® cm, at the appropriate Keplerian 
orbital speed (~ 2.074kms“^). The hydrodynamics of the 
gas was neglected for this test, and the sink particles were 
allowed to move through the stationary gas under gravi¬ 
tational forces only. The benchmark test was run for 130 
orbits of the inner-most particle (corresponding to ~ 26 or¬ 
bits of the outermost particle). The orbits of the particles 
are overlaid on the density distribution and adaptive mesh 
in Figure pTj indicating that the integrator, domain decom¬ 
position, and gas-particle forces are operating satisfactorily. 


7 ACCRETION TESTS 


where Mj is the mass of the sink particle, and Vij and f 
are the distance and direction vector between the cell centre 


The principle tests of our sink particle implementation are 
those that ensure that the accretion of gas occurs at the 
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Figure 10. The paths of three masses in the Pythagorean three-body benchmark test. The panels shows show the paths of the sink 
particles at f = 0 — 10 (left panel), f = 10 — 20 (middle panel), and t = 20 — 30 (right panel). Mass 1 is shown as a green, dashed line, 
mass 2 as a red, solid line, and mass 3 as a blue, dotted line. Solid circles indicate the position of the particles at unit time intervals. 
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Figure 11. The results of the gas-on-sink gravitational force test. 
The paths of the three test sink particles are shown as bold solid 
lines, whilst the density of the gas is shown as a logarithmic colour 
scale. The AMR mesh is shown by thin solid lines. 


correct rate. We therefore conducted three benchmark tests 
of increasing complexity. 


7.1 Bondi accretion 

A point mass M accretes spherically from a gas cloud, where 
far from M the cloud is stationary and has density poo- 
|Bondi| ( |1952[ ) showed that the maximum accretion rate is 
given by 


M = 47rA^:^-^poo (21) 

Cs 

where Cs is the sound speed and A = 1.12 for isothermal gas. 
The corresponding Bondi radius is given by 


Rb 


2GM 


( 22 ) 


For this test we adopted M = 1 Mq and poo = 
10“^® gcc“^, giving an expected accretion rate of 5.9 x 10“^^ 
Mq yr“^and a Bondi radius of 3.75 x 10^^ cm. We ran three 
two-dimensional models with resolutions of Rb/A x of 5, 10, 
and 20. Rapid convergence in the accretion rate with reso¬ 
lution is observed (see Figure 121. The highest resolution 
model has an accretion rate of 5.7 x 10“^^ MQyr”'^, i.e. 
within 2% of the theoretical rate. 


7.2 Collapse of a singular isothermal sphere 


Shu 


(19771 showed that an isothermal sphere with p(r) oc 


1/r^ collapses in such a way that there is a constant mass 
flux through spherical shells. The added level of complexity 
over the Bondi test described above is that the self-gravity 
of the gas is significant. 


We adopt the same test as Federrath et al. (20101, 


with a sphere of radius i? = 5 x 10 ° cm with p = 3.82 x 
10“^® gcm“® with contains 3.02 Mq of gas. The sound speed 
of the gas was taken to be 0.166 kms“^. The expected accre¬ 
tion rate of this model is 1.5 x lO”"^ Mq yr“^. We adopted an 
adaptive, two-dimensional cylindrical mesh for four levels of 
refinement, with the smallest cells of R/Ax = 300 

The model demonstrated excellent agreement with the 
Shu prediction, starting with an accretion rate of ~ 1.5 x 
lO^'^MQyr”^, and only declining when 90% of the origi¬ 


nal mass had been accreted (see Figure 131. At the end of 
the run the local density approached the global floor den¬ 
sity for the simulation (10“^^ gcc“^) and the accretion rate 
approaches the Bondi rate for that density. 


© 2015 RAS, MNRAS 000,[iHTT] 


































































































































































































































































































































10 Tim J. Harries 



Time (ta) 

Figure 12. Results of the Bondi accretion test. The accretion 
rate is plotted against the Bondi timescale {tg = Rg/cs) for the 
Rg/Ax = 5, 10, and 20 models (dotted, dashed and solid lines 
repsectively). The thick solid line corresponds to the expected 
theoretical rate of 5.9 X 10“^^ Mq yr“^. 

7.3 Bondi-Hoyle accretion 



0 5 10 15 20 25 30 

Time (kyr) 


We constructed a 2D test case in which a 1 Mq sink is placed 
at the origin in initially uniform density gas (10“^®gcc“^) 
with molecular weight of 2.33 and a temperature of 10 K. 
The gas initially had a constant velocity with a mach num¬ 
ber of At = 3 parallel to the z-axis, and an inflow condition 
at the upstream boundary and and outflow condition at the 
downstream boundary. These are the same initial conditions 
as the Bondi-Hoyle test case in |Krumholz et aL] ( |200^ , al¬ 
though they ran their simulation in 3D. 

The model extent was 0.78 pc, and three levels of re¬ 
finement were used with the smallest cells corresponding to 
2.3 X 10^^ cm. We ran the model until an approximately 
steady-state of the accretion rate was achieved (Figure [l4|). 
The theoretical Bondi-Hoyle accretion ra te for these con- 
ditions is 1.7 x 10“^^MQyr“^, but the Krumholz et al!] 
(20041 model, and those of Ruffert (19961, found accretion 
rates of close to 2 x 10~^^ Mq yr””^, with considerable tem¬ 
poral variation in the accretion rate on the Bondi-Hoyle 
timescale. Our simulation reached a st eady-state accretion 
rate of 2.4 x 10~^^ Mq yr“^(Figure 151, which is compara¬ 
ble to the peak accretion rate seen in the [Krumholz et ah] 
(20041 simulations. The level of variability is substantially 
lower in our simulation, presumably because the instabilities 
that build up and modify the dynamics near the sink in the 
3D simulations are absent in our 2D models. 


8 CONCLUSIONS 

We have presented a new method for including radiation 
pressure in RHD simulations that incorporates a level of mi¬ 
crophysical detail that is significantly greater than that of 
flux-limited diffusion or hybrid techniques. We have shown 


Figure 13. Results of the Shu accretion test. The lower panel 
shows the accretion rate as a function of time (crosses), with the 
expected theoretical accretion rates predicted by the Shu model 
(solid line) and Bondi accretion at the floor density (dashed line). 
The upper panel shows the growth of the central object as a 
function of time (crosses) along with the expected trend assuming 
a constant Shu accretion rate (solid line). 


that the new method works in both the pure absorption 
and pure scattering regimes, and properly treats anisotropic 
scattering processes. The method comprises a simple addi¬ 
tion to the MC estimators required for radiative and pho¬ 
toionisation equilibrium calculations, and does not therefore 
represent a substantial computational overhead to the MC 
RHD methods described by Haworth & Harries (20121. 

However the MC method as a whole is significantly 
more computationally demanding than the FLD and hybrid 
methods, and we have therefore developed two new methods 
to ameliorate this. The first is the packet splitting method 
detailed in section in which we extend Lucy’s original al¬ 
gorithm to incorporate photon packets with varying energy. 
The second is to distribute the MC photon packet loop over 
many instances of the computational domain, allowing an 
excellent scalability (see section]^. 

The final element needed to compute massive star for¬ 
mation models is a description of the protostar itself, and 
to do this we have included a sink particle algorithm. By 
interpolating on protostellar evolutionary model grids (e.g. 


Hosokawa & Omukai 20091 as a function of mass and ac¬ 


cretion rate will will be able to determine temperatures and 
luminosities for our protostar, and assign it a spectral en¬ 
ergy distribution from appropriate atmosphere models. The 
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Figure 14. Bondi-Hoyle accretion test. The figure shows loga¬ 
rithmically scaled density between 10“^^ gcc“^ (pink) and 10“^® 
gcc“^ (black). The white arrows show velocity vectors, with the 
longest arrows corresponding to speeds of 3kms“^. The black 
semicircular region at the origin signifies the extent of the accre¬ 
tion region. 



Figure 15. Bondi-Hoyle accretion as a function of time. Time is 
given in units of the Bondi-Hoyle timescale, while the accretion 
rate is given in units of 10“^^ Mq yr“^. 


sink particles then become the origin of photon packets for 
the RHD calculation. 

The next stage is to compute massive star formation 
models that incorporate radiation pressure and ionisation 
feedback. Initially these models will be two dimensional, but 
we also conduct three-dimensional calculations in order to 


simulate binary star formation. Of course the algorithms 
detailed here have wider applicability, and the correct treat¬ 
ment of feedback from ionisation, radiation pressure and 
winds is important for cluster-scale calculations, in particu¬ 


lar with reference to gas dispersal from clusters (e.g. Rogers 
fc Pittai^|2013 1. 
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